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Abstract 



We have implemented non-ideal Magneto-Hydrodynamics (MHD) effects in 
the Adaptive Mesh Refinement (AMR) code RAMSES, namely ambipolar diffu- 
sion and Ohmic dissipation, as additional source terms in the ideal MHD equa- 
tions. We describe in details how we have discretized these terms using the 
adaptive Cartesian mesh, and how the time step is diminished with respect to 
the ideal case, in order to perform a stable time integration. We have performed 
a large suite of test runs, featuring the Barenblatt diffusion test, the Ohmic dif- 
fusion test, the C-shock test and the Alfven wave test. For the latter, we have 
performed a careful truncation error analysis to estimate the magnitude of the 
numerical diffusion induced by our Godunov scheme, allowing us to estimate the 
spatial resolution that is required to address non-ideal MHD effects reliably. We 
show that our scheme is second-order accurate, and is therefore ideally suited 
to study non-ideal MHD effects in the context of star formation and molecular 
cloud dynamics. 

Subject headings: stars: formation — ISM: magnetic fields — methods: numerical 
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1. Introduction 

The impact of magnetic fields on various objects in astrophysics is now well established. 
They play a major role on a wide range of scales, from the study of the early universe, the 
stellar and intergalactic medium to the formation and interiors of stars or the accretion 
flows around stellar objects. They are difficult to study both from an observational and 
a theoretical (and numerical) point of view. Several implementations of ideal MHD have 
been performed since the last decade (Fromang et al. (2006), Stone & Norman (1992), 
Machida et al. (2005) among others), and numerical issues concerning the divergence free 
condition have now been resolved. However, ideal magnetohydrodynamics (MHD) is in 
many circumstances a poor approximation and non-ideal MHD effects need to be thoroughly 
considered. 

Ambipolar diffusion is expected to play a major role in star formation (Mestel & Spitzer 
(1956)), at the scale of molecular clouds by enabling the collapse of otherwise magnetically 
supported clouds (Basu & Ciolek (2004)) and at the scale of the first Larson's core with 
the formation of a centrifugally supported disk and the well-known fragmentation crisis 
(Hennebelle & Teyssier (2008)). Ambipolar diffusion is also important in protoplanetary 
disks as they are in general only partially ionised. The microscopic and entropic heating 
resulting from the drift and collision between neutral and charged species is another very 
important and relatively unknown aspect which is crucial as soon as cooling or heating of 
the gas (thus radiative transfert) is taken into account (in contrast it is not relevant when 
using a barotropic equation of state). 

Magnetic resistivity effects range from prohibiting long-term MHD turbulence in 
molecular clouds (Basu & Dapp (2010)) to preventing the magnetic braking catastrophe 
on small scales (Dapp & Basu (2010)). Its importance is also crucial in order to study 
disk formation around protostellar objects (Krasnopolsky et al. (2010)) and the physics of 



-4- 



binary formation and brown dwarfs. 

Thus, it appears necessary to introduce the ambipolar and Ohmic diffusion in a 3D 
MHD code. Before exploring the astrophysical impact of such a study, however, the 
accuracy of the treatment of the complete MHD set of equations must be unambiguously 
assessed. This is the very aim of the present paper, in which we describe a prescription 
to incorporate ambipolar and Ohmic diffusion in the multi-dimensional MHD AMR code 
RAMSES (Teyssier (2002)), extending the ideal MHD version presented in Teyssier et al. 
(2006) and Fromang et al. (2006). 

Several numerical treatments have been derived from ideal MHD models, following 
different aims and thus using different methods. The first attempt to implement ambipolar 
diffusion in a code was made by Black & Scott (1982) using an iterative approximation in 
an implicit first-order code. Toth (1994) used a semi-explicit method in a two-dimensional 
code to investigate instabilities in C-schocks. Mac Low et al. (1995) presented a widely 
used explicit method (Choi et al. (2009), Mellon & Li (2009), Li et al. (2011)) to implement 
single- fluid ambipolar diffusion in the strong coupling limit, and then developed a two-fluid 
model in order to capture shock instabilities. Tilley & Balsara (2008) and more recently 
Tilley & Balsara (2011) presented a semi-implicit scheme for solving two-fluid ambipolar 
diffusion, arguing that the single fluid approximation does not carry the full set of 
MHD waves that can propagate in a poorly ionized system. Multi-fluid approaches 
including ambipolar diffusion and Ohmic diffusion have been suggested by Falle (2003), or 
O'Sullivan & Downes (2006) and then investigated by e.g. Kunz & Mouschovias (2009). 
Recently, Li et al. (2011) used the single-fluid approach including more realistic resistivities 
based on a multi-fluid approach for ambipolar diffusion, Ohmic diffusion and Hall effect 
in two-dimensional (axi-symmetric) calculations. Another approach has been used by 
Machida et al. (2006) was to describe both ambipolar diffusion and Ohmic diffusion in one 
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single Laplace operator r]AB, with 77 taking into account every diffusive process at stake. 
These numerous studies have also given rise to several numerical tests, a number of which 
we will either perform directly or slightly modify to assess the accuracy of our treatment. 

Our current study focusses on the one-fluid approximation (Shu et al. (1987)), as in 
previous calculations by Mac Low et al. (1995) and Duffin & Pudritz (2008). We used a 
direct explicit method to implement non-ideal MHD terms in both the induction and energy 
equations (Mac Low et al. (1995)) in an AMR framework. We did not choose to account for 
non-ideal effects by adding ambipolar diffusion and Ohmic dissipation in a single Laplace 
operator as Machida et al. (2006). Instead we kept the full expressions and proceeded 
separately for each non-ideal effect. 

The paper is organized as follows. In § 2, we first derive the equations for ambipolar 
diffusion in the single fluid approximation. We then describe the various tests we have 
performed, first without the hydrodynamics and then in a complete MHD situation, 
exploring in particular the propagation of Alfven waves. Comparisons with existing 
analytical or benchmark solutions are presented in details, demonstrating the validity and 
the accuracy of our scheme. § 3 addresses the case of Ohmic diffusion, following the same 
procedure as for ambipolar diffusion, while § 4 is devoted to the conclusion. 

2. Ambipolar diffusion 
2.1. Equations 

When the ions pressure and momentum are negligible compared to those of neutral 
species (as is the case for example in molecular clouds), the Lorentz force exerted on the 
ions is in equilibrium with the drag force exerted by the neutrals, which corresponds to 
a situation of strong coupling between the neutral fluid and the field lines. In such a 
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situation, the plasma can be adequately described by a single fluid (Shu et al. (1987) and 
Choi et al. (2009)) of mass density p, and neutrals and ions mass densities p n ~ p and 
Pi <C p n respectively. Interestingly, in the case of the one-fluid approximation the results 
can be directly compared with the ones obtained with ideal MHD giving clear insights 
about MHD wave propagation properties in the non-ideal case (Balsara (1996)). The present 
study is devoted to the technical resolution of the resistive MHD equations and we will also 
ignore gravity (and thus the Poisson equation). The MHD equations are given by the usual 
continuity, momentum, energy and induction equations, completed by the magnetic field 
divergence-free condition (in rational units, B rat = B CffS /v^4vr): 



f + V.( P v) = (1) 

dv 

p— + p(v.V)v + VP-F L = (2) 

+ V.{(E tot + P tot )v - (v.B)B - E AD x B) = (3) 
<9B 

--Vx(vxB)-VxE iO = (4) 

V-B = 0. (5) 



E tot denotes the total energy 



E tot = pe+^pv 2 + l -B\ (6) 



where e is the specific internal energy. 
P to t is the total pressure 



P to t = {l-l)pe+ l -B\ (7) 



where 7 is the adiabatic index. 
F L denotes the Lorentz force 

F L = (VxB)x B, 
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with 



1 



(9) 



v ; - v n = 



F L . 



lADPiP 



The ambipolar electromagnetic force (EMF) is given by 



Ead = (Vi — v„) x B 



1 



(10) 



F L x B, 



lADPiP 



where v$ and v n denote respectively the ions and neutrals velocities, and ^ad is the drift 
coefficient between ions and neutrals due to ambipolar diffusion. The last equality in 
Equation (10) illustrates the balance between magnetic and drag forces in the ion fluid, 
while Equation (3) is accounting for ambipolar heating of neutrals by ions (Shu (1992)). In 
order to write both Equation (3) and Equation (10) we need to assume that the velocity 
drift between ions and neutrals and the one between electrons and neutrals are not too 
different (by a factor ~ Therefore, as long as the Hall effect is negligible, these 

equations remain valid (see Pinto et al. (2008) for a more detailed study). 

Equation (3) (the conservation of energy: + V ■ Energy = 0) is equivalent to 



where we can see that the neutrals-ions friction term heats up the gas and increases the 
entropy. 



PT-JT 



ds ||(VxB)xB 



(11) 



dt lADPPi 



2.2. Computing the ambipolar diffusion terms 



In this section, we describe the numerical implementation of the previous equations, 
focusing on the ambipolar diffusion terms in the energy and induction equations. 
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2.2.1. The ambipolar EMF 

The ambipolar term in the induction equation can be considered as an additional 
electromotive force (EMF). To update the magnetic field the values of the EMF have to be 
defined as time and space averages along cell edges (Teyssier et al. (2006)). For instance, 
the EMF in the z direction is defined at x^i , y^i, z\. (and the same for the other directions 

2 •* 2 

with circular permutations) where i,j and k are the cell indices in the x, y, z directions 
respectively. 

We focus now on the EMF in the z direction, defined at x { _i , j/ 7 -_i, Zk and explain in 

2 •'2 

details how it is computed. writes 

E AD = — - — [(V x B) x B] x B = v d x B (12) 

lADPiP 

with the drift velocity vh = v; — v n = — - — F L . We therefore have to evaluate — - — , F L 

J U 1 " lADPiP lADPiP 1 ^ 

and B at x,_i , , 2%, and then calculate 

2 J 2 

The RAMSES code is based on the Constrained Transport scheme for the magnetic 
field evolution (Teyssier et al. (2006)), for which the components of the magnetic field are 
defined at the center of cell faces: if Xi,yj,Zk are the coordinates of a cell center, B x is 
defined at position x { _^y^z k , B y at Xi,yj_i,z k and B z at Xi,yj,z k _i (see Figure 1). Each 
magnetic field component is computed using the finite-surface approximation, which reads 
for the x component: 

\ i i r y i+h [ z >+h 

B x;i-i,j,k) = A y A z J J B x (x i _i,y,z)dydz, (14) 

while other components are defined by circular permutations. 

We also need to define the drift velocity v<j at x^i, yj_i,Zk using the Lorentz force Fl, 
the density p and the ions density pi. The two latter quantities are cell-centered quantities 
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'in contrast to the magnetic field): 



Ax Ay Az 

and 



111 [ x *+h f y *+h 
(Pi,j,k) =^—^—^— / / p(x,y,z)dxdydz, (15) 



x. i Jy. i Jz. i 



ill r>+i f y i+h f z '+i 

(Pions;t,j,fe) = ~^ /±y/\ z J J J pi(x,y,z)dxdydz, (16) 



z. i 

2 



We then define the edge-centered density (and the ions edge-centered density) as an 
arithmetic average of surrounding cells (see Figure 1, right panel): 

Pi-yik =-j\Pij,k + Pi,j-i,k + Pi-i,j,k + Pi-x,j-i,k\- (17) 



This definition is adapted for the components of the magnetic field to compute them 
at x i _i,yj_i,Zk (see Figure 2): 



5, 



2 '■> 2 > K 



2 L 
1 



B 



+ 5 



x;i—-x,j—l,k 



(18) 
(19) 



Given the Lorentz force at , y , zj. (see Sections 2.2.2 and 2.2.3), it is possible 
to compute the first component (z direction, with unit vector ez) of the ambipolar EMF: 
Ead ' e z — ._ 1 k , while the two other components are obtained through circular 

permutations. 

These ambipolar EMFs are then added to the ideal MHD EMFs calculated with a 2D 
Riemann solver, as described in Teyssier et al. (2006) and Fromang et al. (2006). 
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i-l,j,k 
+ 



i-lj-l.k 
+ 



i,j,k 

Jx;i-l/2,j,k 



>y;i.j-i/2,k 



MR,. 



i,j-l,k 
+ 



n-ij,k 
+ 

1 


1,k 
+ 

O 

^1-1/2 j-l/Z,k 


^-l,j-l,k 





Fig. 1. — Left: coordinates of the center of neighbouring cells, with the natural places 
where the magnetic field and the EMF are defined. Right: computation of the density at 



Zk by averaging over neighbouring cells. 



i-l,j,k 
+ 


Bx^i-l/2,j,k '/J>k 
1 

1 

.1 B x ;i-l/2,j-l/2,k 


i-l,j-l,k 
+ 


— >T 

1 

\i,j-l,k 
B * + 



VX 



i-l,j,k 
+ 

By;i-l,j-l/2,k —J 

r" 


i, 

H 


j.k 

By ; |,j-l/2,k 


i-l.j-l.k 
+ 


1-1/2 ,)-l/2,k 

i,j-l,k 
+ 



VX 



Fig. 2. — Computation of the magnetic field B x and B y at x i _i,yj_i, Zk- 
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2.2.2. The Lorentz force as the product of the current and the magnetic field 

We now focus on the computation of the Lorentz force Fl at z k . The first 

way to calculate this term is to explicitly compute the magnetic field components and the 
current J = V x B at x t _i , y -_i , z k , as: 

2 J 2 

F L = J x B. (20) 

OB f) R 

Jz = -Qjr — ^f is naturally defined at , Vj_i, z k and J x and J y are naturally defined 
respectively at Xi, z k-\ an d x i-±> Vji z k-\- I* 1 order to have all three components of J 
at the location of the EMF (x^i , j/ 7 -_i , Zk) we need to use the magnetic field components 

2 •* 2 

at specific positions, as follows: 

J* = A* '(/',. ~ ,., ,,/) - ^z-\B mi _ k ._ hk+k - B y ._y_ hk _ h ) (21) 
Jy = ^\B x;l _ hj _ hkH - B x;i _ k!j _ hk _0 - Ai- 1 ^!, - B^-fr) (22) 
J, = A.r 1 - ,., ,.,)) - Ay" 1 (/^,,, - ,,) • (23) 

As above, we have to express each of these terms through arithmetic averages of the 
naturally defined components of the magnetic field: 

;.,./,■ and B y ..._ hk (24) 

B z;i,j-±,k aIld B z;i-\,j,k (25) 

^.i-ij-i.fc-i and B^i^i^i. (26) 

and are naturally defined at x t _\.yA.Zi. and Xi,y,_i,z k respectively, and thus 
the two terms in Equation (24) need not to be computed again. 

The terms in Equation (25) and Equation (26) are obtained thanks to averages on four 
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components (see Figure 3): 



^z;i-i,j,k ~ 4 [Bz;i,j,k+i + B z\i,j,h-\ + & 'z;i-l,j,k+~ ^ ^ z;i-l,j,k-\\ 

= ^ [Bz;i,j,k + -Sz;i-l,i,fe] (27) 

B z . itj _i k = - [B z . iJtk + Bj-ij-iJ (28) 

S a , i _i ij _i >A ._i = -[B^ij^ + B^ij^ + B x;i _i J)fe _ 1 + S a ,. i _i J _ life _ 1 ] (29) 



The Lorentz force for the other components of the EMF (x and y directions) are 
obtained through circular permutations. 

2.2.3. The Lorentz force as the divergence of a flux 

Another way to compute the Lorentz force is to express it as the divergence of a well 
chosen flux: 

F L = (V ■ .FOei, (31) 

and 

J~\ = BiBj&j <5jj6 ma g6j, (32) 
with i,j E [x, y, z] and e mag = \{B 2 X + B 2 y + B 2 Z ). 

Let us focus on the x component of the Lorentz force for the EMF in the z direction. 
It reads: 

F L • e x = d x (B% - e mag ) + d y (B x B y ) + d z {B x B z ). (33) 
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In order to compute the Lorentz force at , y , Zk, we compute each term at 

2 -'2 

specific positions: 





B 2 , i J 








^;i-lj-2,fe)] 


dy{B x B y ) =—[B X f_y 









(34) 
(35) 

(36) 

We then only need to compute each component of the magnetic field at x i _i,yj_i,z k 
in order to get the EMF in the z direction. 

As explained in the previous paragraph, an average over well chosen (where the 
magnetic field is naturally defined) surrounding cells is used (see Figures 2 and 4): 



B, 



B x , 




B z 




By 




hi- 


2' K 2 



\_Bx\i,j,k ~\~ ■Bx;ij—'\.,k\ 
[B Z ;i,j,k + B z -ij_i jk \ 
\By;i,j,k "t - By-i— l,j,k\ 



[B X ;i-\,j,k + ^x;i-i,3-l,fe 



B x; i-i t j t k-l + ^x;i~2 ,j-l,fe-l] 
1. 



(37) 
(38) 
(39) 

(40) 
(41) 
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and 







2 I* 1 




y- 


y 




y- 


y 



J B y;i,j-±k + By^ji 



2> fcJ 



4 
1 

2' 

1 

2' 

1 

2' 

1 



7j( B z;i,j,k~± + B z;i,j-l,k-i) 



z:i- 


-l,i,fc-i 


+ B z;i-l,j-l,k-± 


z;i,_ 


i,k+± + 


B z;i,j-l,k+±) 


z;i- 




+ B z;i-l,j-l,k+± 



[-^z;i,j,fc "I - B z;i,j—1, 



(42) 
(43) 



(44) 



Again and as highlighted previously, in order to get the two other components of the 
EMF one only needs to perform circular permutations. 

These two methods (described in Sections 2.2.2 and 2.2.3) to compute the Lorentz force 
are implemented in RAMSES and show similar performances. This method might work better 
under certain conditions, for a particular setup of the magnetic field lines. Nonetheless, 
when counting the number of floating point operations, the computer performs using this 
method 4911 more additions and 8047 more multiplications for a given cell than with the 
previously described method. 

2.2.4- Contribution of ambipolar diffusion to the energy flux 

The ambipolar energy flux (see Equation 45) has to be evaluated on each face of the 
cell, that is to say at locations (x i± i,yj,z k ), (x iy y +i , z k ) and (a?j, yj, z k± i). Again, as in 

2 ** 2 2 
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3 z;i-l,j,k+l/2 




B X ;i-l/2,j-l,k 
( 


B X !!-l/2,j,k 
Bx;i-l/2,j-l/2,k-l/2 


R® 

D x;i-l/2,j-l,k-l 


B® 

D x;i-l/2,j,k-l 



Fig. 3. — Left: B z;i _i - k as an average over surrounding cells. Right: B xi _i -_i k _i. 
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] 2.2.1, the needed components are obtained thanks to averages over neighbouring cells 
(averages which are not detailed here). 

jr AD = _e ad x B = — ((J x B) x B) x B. (45) 

lADPiP 



2.2.5. Computation of the time step in presence of ambipolar diffusion 

The ambipolar diffusion timescale can be estimated through the drift velocity of ions. 
Recalling Equation (10) we get: 

||vdrift|| oc ||F L || 

(46) 

lADPi^AD 

where Lad is a characteristic length for ambipolar diffusion, which can be estimated as 

^ad = V ||^^ ■ We then have the timescale: 

Lad ^ADpiL\ D . . 

tad = 7i rr = o • I 47 ) 

||v d rift|| V A 

L 2 v 2 

Written as a diffusion, tad — with the ambipolar diffusion coefficient D — — — 



D r *YADPi ' 

where va = ^= is the Alfven speed and {^adPi)~ 1 is the characteristic collision time 
between ions and neutrals. A Von Neumann analysis for the diffusion part of the equation 
can be performed for the scheme used: 

<9B 

— - V x E AD = 0. (48) 



It can be differenced (in one dimension): 

rm+l r>n 



Using m = e n e ikjh , 



D/\t, 

e = 1 + 2 T {cos{kh) - 1). (50) 

Ax 
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Equation (50) shows that the scheme is stable according to Von Neumann stability analysis 
provided the coefficient is lower than 0.5: 

lAx 2 lAx 2 -f A DPiP / C1N 
\e\<l*At< — = - B2 . (51) 

For the three dimensional case, this time-step constraint is more stringent than for the one 
dimensional case presented above. 

Therefore, the time step used to update the solution is computed by taking the 
minimum of the usual MHD Courant condition (Fromang et al. (2006)) and the ambipolar 
timestep defined by 

t AD = 0.1 x min(^^ Ax 2 ), (52) 

V A 

where the minimum is taken over all the cells of the computational grid. The coefficient 
0.1 < | is taken to achieve better convergence. This choice is based on the various tests 
performed, and might not be suited to all other problems. As can be noted in equation (52), 
the time-step scales as Ax 2 . Even though this is very demanding in terms of numerical 
resources as the grid becomes more and more refined, there are means to speed-up the 
calculations, as explained in the following paragraph. 

The ambipolar time step is proportional to pi (see Equation (52)), which is assumed to 
be proportional to p k : pi = C^fp (see Elmegreen (1979)). Both the factor (C) and the power 
law (pi) are very dependent on the microphysics and the geometry of the grains. This 
assumption is thus made for the sake of simplicity, but might not always be valid. In some 
cases, for example in star formation simulations, the time-step can become unphysically 
small in very diffuse regions where the ionisation approximated as above (Equation (52)) 
is very small. Following Nakamura & Li (2008), we use a threshold in order to limit the 
time-step when needed. On the other hand, in very dense parts where the grid is fairly 
refined (where Ax is small), the dependence of pi and ^ad with the density prevent the 
time-step from becoming too becoming too small. This situation has to be studied for each 
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different physical problem and can't be assumed once and for all. We will address this issue 
in the case of star-formation in a forthcoming paper. 

2.3. The AMR scheme 

The AMR algorithm used in RAMSES is described in Teyssier (2002), and its extension 
to MHD is first described in Teyssier et al. (2006) and then in Fromang et al. (2006). We 
briefly recall the main features here. It is a tree-based AMR code whose data structure is a 
'Fully Threaded Tree" (Khokhlov (1998)). The grid is divided into "octs" which are groups 
of 8 cells with the same parent cell. The first level of refinement (/ = 1) corresponds to the 
unit cube, which defines the computational domain. The grid is recursively refined from 
the / = 1 to the minimum level of refinement l m in, in order to build the base Cartesian grid. 
Adaptive refinement then proceeds from this coarse grid up to the user-defined maximum 
level of refinement l max . When l max = l min the computational grid is a traditional Cartesian 
grid. Issues arise when refined cells are created, in the case where l max > l min . Concerning 
the non-ideal MHD, the EMF and energy fluxes are simply added to the existing ideal 
MHD EMF and fluxes. As a consequence, there are no more complications in refining and 
derefining cells than in the ideal MHD case. 

2.3.1. Divergence-free prolongation operator: refining cells 

The "prolongation operator" is the creation of a new "oct" of 8 cells when a cell is 
newly refined. Cell-centered variables and magnetic field components are needed for each 
refined cell. This is usually done using a conservative interpolation of the variables, yet in 
the case of magnetic fields, the divergence-free constraint has to be fulfilled by each of the 
new cells which makes things more complicated in details. A critical step has been solved 
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by Balsara (2001) and Toth & Roe (2002) in the constrained transport framework. The 
idea developed in those articles is to use slope limiters to interpolate the magnetic field 
components in each parent face conserving the flux, and then to perform a three dimensional 
(which is divergence-free inside the cell volume) reconstruction in order to compute the new 
magnetic field components for each children faces. The same slope limiters as the ones used 
in the Godunov scheme for the hydrodynamics are used in this step. 

2.3.2. Magnetic flux corrections: derefining cells 

The " Restriction Operator" is, in the multigrid terminology, the operation of derefining 
a split cell. The divergence-free constraint still needs to be satisfied, so that the magnetic 
field components in the coarse faces are simply the arithmetic averages of the four fine faces 
values. This is the parallel in MHD of the "flux correction step" for the Euler system. 

2.3.3. EMF corrections 

This is specific to the induction equation: for a coarse face adjacent to a refined face, 
the coarse EMF in the conservative update of the solution needs to be replaced by the 
arithmetic average of the two fine EMF vectors. This is mandatory to guarantee that the 
magnetic field remains divergence- free even at coarse/fine boundaries. 

2.4. Tests for the ambipolar diffusion 

2.4-1- The Barenblatt diffusion test 

In this section, we first test the accuracy of the calculation of the ambipolar term alone. 
For sake of simplicity, we assume that the magnetic field has the form B y (x, z), with B x = 
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and B z = 0, that all the velocities remain equal to zero and that density and thermal 
pressure are constant. The induction equation takes the form of a diffusion equation: 

% d x | d / Bl 
dt dx\^ADPipdx V J dz\<y AD pipdz 

which can also be written in compact form: 

dB„ _ / v 2 



^ = J.(_S_» (fl>) ) + |.(_S_a w ), (53 ) 



V.(-^-vV). (54) 



This is a non- linear diffusion equation, since the diffusion coefficient, tjad 



IAD Pi 1 

depends non-linearly on the magnetic field. Here, va = B y j yfp denotes the y-component 
of the Alfven velocity. The solution of this problem with a Dirac pulse as initial condition 
(known as the Barenblatt-Pattle solutions) has been derived by Grundy & McLaughlin 
(1982) (See Appendix A for more details about the analytical solution). 

The initial states in one and two dimensions are respectively: 

1 if X ^ 0.9 Aji ere ; = 3 

B y0 = < (55) 

elsewhere 

1 if \f{^ •^center) ^center) ^ 0.9 Al; cve j = 3 

Byo = \ (56) 
elsewhere. 

with Axi eve i = z being the cell size at the lowest level of refinement used (in this case: 3). 
This ensures that the initial perturbation is the same for the case of an AMR grid and a 
fully refined grid. 

We performed the test both with a uniform grid and using the AMR with the same 
maximum level. The level of refinement refers to the number of cells used: 2 nd cells are 
used for the level of refinement N in a D-dimensional calculation. As seen in the figures, 
the agreement between the numerical and the analytical curves is excellent, a few tenths of 
a percent on average. The results obtained on an AMR grid (with levels varying from 3 to 
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7, corresponding to a mesh size Ax = 0.5 3 and Ax = 0.5 7 ) are almost as good as the ones 
obtained on a regular grid corresponding to the highest level of refinement (level 7, with 
a cell size Ax = 0.5 7 and 128 cells): the maximum relative error is less than one percent 
except where the magnetic field equals zero. The difference between AMR and uniform grid 
is less than 2.10~ 4 for values of magnetic field of about 0.01. 

The results for B y (x) are shown on Figure 5, where we have taken jad — 1> Pi — 1 and 
p = 1 and on Figure 6 for B y (x, z). 

The grid is refined if the gradient of magnetic field is greater than 0.1 (this insures for 
this test that the error on the AMR grid and on the regular grid are about the same). We 
also checked that the same accuracy is obtained for any orientation of the magnetic field 
(using Bx or Bz instead of By). Figure 7 represents the evolution of the error calculated 
as e = y^E^Ii ^ umCTical ^^ nal iM££j! as function of the mesh size (N being the number of 
cells for each level). 

In terms of computational time for this particular test using the refinement strategy 
described above, the time is about the same for a regular grid at level 7 as for a grid going 
from levels 5 to 7, but there is a gain of about 40% in the number of cells. One level further 
(regular grid at level 8 or AMR grid going from levels 5 to 8) the computation is 30% faster 
in the AMR case and needs 55% less cells. For a regular grid at level 9 or an AMR grid 
going from levels 5 to 9 the calculation is 60% faster with 70% less cells needed. 

2.4.2. The C- shock test 

Following Duffin & Pudritz (2008) and Mac Low et al. (1995), we have tested our 
new scheme for the case of both isothermal and non-isothermal oblique C-shock including 
ambipolar heating as given by Equation (3). We start from a steep function as initial state 
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Fig. 5. — Barenblatt diffusion test for ambipolar diffusion at t = 200 with B y being a function 
of x only. The left panel is a snapshot of the AMR run with levels from 3 to 7. The right 
panel corresponds to a fully refined Cartesian grid up to level 7. 
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Fig. 6. — Barenblatt diffusion test in 2D with B y depending on x and z. Here, the calculation 
was performed on an AMR grid from level 2 to 6. The left snapshot is a 2D contour plot at 
t = 200: the symmetry of the solution is preserved. The right snapshot (same legend as in 
Figure 5) is a ID cut across the maximum at t = 200. 
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for the different variables, whose values are the ones taken at infinity ahead of and behind 
the shock. Our calculation takes place in the frame of the shock. The post- and pre-shock 
values are displayed in Table 1. The angle between the shock normal and the magnetic field 
is set to 45°. 

For this test, we set jad = 75, pi = 1. The sonic Mach number is M, — 10 and the 
Alfven Mach number is M.a — 1-8. Outflow boundary conditions are used in the simulation. 
After a short transient phase, the shock becomes stationary. 

The isothermal shock is modeled through P n = p n c 2 s with c s = 0.5 the sound speed, and 
without solving the energy Equation (3). Results are shown in Figure 8 and compared to 
the semi-analytical solution described in Mac Low et al. (1995) (see Appendix B for more 
details). 

For the non-isothermal case the energy Equation (3) is solved assuming a perfect gas 
with an adiabatic index 7 = | and without any additional cooling. The semi-analytical set of 
equations to be solved is derived from Duffin & Pudritz (2008), where we assume a constant 
ion density (see Appendix C for more details). The steady-state is not very different from 
the isothermal case, except for the pressure. The results for the non-isothermal case are 
shown Figure 9. Our results are significantly different from Duffin & Pudritz (2008) in the 
pressure across the shock. This is explained by the additional heating term (and an artificial 
cooling term necessary for the equations to converge) in their set of equations. Therefore 



Variable 


P 


v x 


Vy 


B x 


By 


P 


Pre-shock value 


0.5 


5 





V2 


V2 


0.125 


Post-shock value (isothermal) 


1.0727 


2.3305 


1.3953 


V2 


3.8809 


0.2681 


Post-shock value (non-isothermal) 


0.9880 


2.5303 


1.1415 


V2 


3.4327 


1.4075 



Table 1: Initial conditions used for the oblique C-shock test, as described in Section 2.4.2. 
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Fig. 8. — Isothermal oblique shock with ambipolar diffusion. Lines and symbols are the 
same as in Figure 5. The levels of refinement vary from 5 to 7. 
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the equations tested are not exactly the same and thus neither are the semi-analytical 
solution nor the results. 

In astrophysical simulations solving non-isothermal ambipolar diffusion only makes 
sense if cooling or heating of the gas is properly taken into account, i.e. if radiative transfer 
is solved. Otherwise, the set of MHD Equations ((1), (2), (4) and (5)) is closed by an 
equation of state (a barotropic one in most cases). 

We also checked that the results are similar for any orientation of the initial magnetic 
field and velocity field. Using AMR gives results almost as good as with a regular grid 
corresponding to the highest level of refinement (not displayed here for conciseness). 

The grid is refined if the gradient of magnetic field, pressure, density or velocity is 
greater than 0.1 (this insures for this test that the error on the AMR grid and on the 
regular grid are about the same). 

2.4-3. The Alfven wave test 

Studying the decay of Alfven waves in an ionized plasma provides a stringent test of 
the coupling between the flow and the magnetic field due to ambipolar diffusion. Following 
Choi et al. (2009), we have examined the behaviour of propagating and standing Alfven 
waves in such a plasma. We closely follow the prescription and the notations defined by 
Lesaffre & Balbus (2007) for the study of Alfven waves in a plasma with Ohmic diffusion, 
and adapt them to the ambipolar diffusion case. Here we derive exact solutions for torsional 
Alfven waves in a non-isothermal plasma with ambipolar diffusion. 
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Fig. 9. — Non-isothermal oblique shock with ambipolar diffusion. Lines and symbols are the 
same as in Figure 5. The levels of refinement vary from 5 to 7. 
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The unperturbed state without Alfven waves is defined as: 



Po = l,Pon = 1,/Oot = 1 (57) 
K,x = 0,K), = 0,^ = 
Box = 0, -Boy = 0, B 0z = 1. 

We seek for perturbed solutions of the form 

u = <5uexp (st + ikz) (58) 
and B = B z + b = B z + <5bexp (st + ikz), (59) 

where 5h = 5b x x + 5b y y and 5u = 5u x x + 5u y y. s is the wave angular frequency and k 
the wave number. For a perturbation wavelength A along the z direction, the wave vector k 
is set to k = 2tt/X. For such solutions, the mass density remains constant along the wave 
trajectory (p = p ). 

Following Lesaffre & Balbus (2007), we restrict ourselves to MHD flows satisfying 
V(P + \B 2 ) = 0, so that the momentum equation reads: 

8 t u = —8 z b (60) 
Po 

and the induction equation simplifies to: 

d t h = B d z u + B ° d 2 z h. (61) 

lADPiOpO 

Combining Equation (60) and Equation (61) gives a quadratic dispersion relation: 

s 2 + k 2 r] AD s + k 2 v\ = 0, (62) 

where the ambipolar diffusion coefficient is defined by t\ad = v 2 a /^adPio and the Alfven 
velocity by v\ = B /^p^. This equation is similar to the dispersion relation obtained 
by Balsara (1996), but we have derived it for the more general adiabatic, non-isothermal 



-29 - 



case, and also for any amplitude in \5h\, provided that V(P + \B ) = 0. If we restrict 



2' 

ourselves to circularly polarized waves with e.g., 5b y = i5b x then V(\B 2 ) = will also 
ensure V(P) = as we now demonstrate. 

It is clear from Equation (62) that Alfven waves propagate (s, ^ 0, with Sj the 
imaginary part of s) only for va > kr]AD/2- The solutions of Equation (62) are given by: 



iz^M-^^pj 2 . (63) 

In the numerical tests that follow, we restrict ourselves to A = 1 and equal to the 
box size, so that k = 2ir. We will explore first a value ^ad — 80, yielding a diffusion 
coefficient tjad — 1-25 x 10 -2 , and resulting in a moderate damping with imaginary part 
Si = ±6.2783387 and real part s r = —0.2467401. We then consider the case 7ad = 30 
{Vad = 0.0333), resulting in a stronger damping with Sj = ±6.2486389 and s r = —0.6579736. 

Estimating numerical diffusion In order to estimate the quality of our numerical 
solution we need to compute the leading order error term in the ideal MHD scheme. This 
is done usually using the Modified Equation approach where a Taylor expansion of the 
numerical solution is performed. We restrict our analysis to the propagation of Alfven 
waves since the Modified Equation is much simpler to handle in this case. We use the 
characteristic variable q; ± = u =)= b/y/p^, so that the system describing the propagation of 
Alfven waves becomes 

9 t a ± ± v A d z o^ = 0. (64) 

We consider here only the right-propagating wave, dropping the superscript ±. The 
conservative update writes 

n+1 _ a n a I - a J 
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Since the Riemann solver accounts for Alfven waves, the interface flux is given by the 
upwind value, solution of the predictor step. 

This entirely defines our second-order accurate numerical solution. We assumed here that 
the time-step is much smaller than the Courant time step, so that VAAt/Az <C 1. Taylor 
expanding the solution and its spatial derivative to the first non vanishing order in respect 
to a™ leads to the following Modified Equation with a second-order leading error term 

d t a + v A d z a ~ d 3 z a ~ ri num d 2 z a. (67) 

The right-hand-side represents a third-order derivative of the solution, usually 
interpreted as a dispersive term. We nevertheless restrict ourselves to the test case studied 
in this paper, namely a sinusoidal wave of period equal to the box size L, and approximate 
the leading-order term as a diffusive term with numerical diffusion coefficient, namely: 

ZkvaAz 2 

Vnum = ' ^ ^ 

From this analysis, we can estimate the amplitude of the diffusion due to the hyperbolic 
solver that needs to be added to the physical (whether ambipolar or Ohmic) diffusion to 
interpret the numerical solution. We also conclude that the leading order term coming from 
the ideal MHD solver scales as Ax 2 . This sets the physical range of ambipolar and Ohmic 
diffusion one can expect to explore for a given mesh resolution. For a mesh of 16 3 cells 
the numerical diffusivity is six times smaller than the ambipolar diffusion with ^ad — 80: 
Vnum = 0.002 and t]ad = 0.0125. This is a good test case in order to assess the accuracy 
of the correction: the dominant term is still coming from the physics, but the numerical 
contribution is not negligible. 

For Alfven standing waves, the same study can be done. Considering two waves: a + 
and a~ , one propagating to the right and the other to the left. The system describing the 
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standing Alfven waves is 

d t a + + VAd z a + + d t a~ — VAd z a~ = 0. (69) 
The interface flux are given by the upwind value for a + 

<| = < + (d z a)t ^ (70) 
= <-i + (71) 

and the downwind value for a~ 

<jf = "J+i - IT (72) 

1 2 Z 

where we then express each term (values and spatial derivatives) in terms of af, using a 
third order Taylor expansion in Az. 

We then obtain for the two propagating waves: 

d t a + + v A d z a + ~ + V -4fdl* + - ^^« + (74) 

O - o _ V A AZ 2 3 _ V^A* 3 4 _ 

-v^a ~ 12 - 48~ ' ^ 

Combining those two equations in order to obtain Equation (69) leads to the solution: 

d t a + + v A d z a + + d t a~ — v A d z oT ~ — — d A z a~ . (76) 

Again, we interpret this fourth order term as a diffusive term with numerical diffusion 
coefficient: 

_ 27TV A AZ 3 
Vnum — 24L 2 ' 

These two expressions for numerical diffusion (oc Az 2 for propagating waves, and 
oc Az 3 for standing waves) are representative of the real diffusion, as confirmed by the 
study of the evolution of the error (Figure 10). 
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To take into account this numerical diffusivity, we solve again the equations of induction 
(Equation 61) and momentum (Equation 60) for a dispersion equation with an additional 
(numerical) diffusion: 

d t u = —d z b + r] num d 2 z u (78) 
Po 

d t b = B d z u + B ° d 2 z h + 7] num d 2 z h (79) 



'jADPiOPO 



yield 



s = _k\y AD + 2r ]num ) ± / _ / JfruoV (go) 



As we have restricted the numerical effect to a diffusion, there is no contribution to the 
imaginary part of the pulsation, as can be seen in Equation (80). 



The propagating Alfven waves test We start the simulation with an initial perturbed 
state with B lx = Ke(6b x e ikx ), Sb x = 1, B ly = 1le(i5b x e ikx ) , and v lnx = TZe(^B lx ) 
and V\ ny = TZe(^^-Biy), where IZe denotes the real part of a complex number. For the 
propagating wave test, we have chosen our initial conditions so that Si > 0. 

The internal energy equation (see Shu (1992)) can be written as: 
dpe . nr7 ((VxB)xB) 2 

-£- + V.(/oev) = -PV.v + ^ '-. 81 

Ot lADpiP 

In the case of perfect gases, we have P = (7 — l)pe. Since Alfven waves are transverse 
waves, V.v = and V.(pev) = 0. The energy equation thus reduces to 

BP 'v — 1 

V = - — -((V x B) x B) 2 (82) 

dt -lADP l p yy 1 ' ^ J 

. This last equation, combined with our choice 5b y = i5b x , gives V(P) = 0. 
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Using Equation (82), the time evolution of the pressure writes 

P = ^ + (7-l)^^(e 2 ^-l). (83) 

Figure 11 shows profiles of Bi x , V\ nx , Bi y , Vx ny , p and P along the z direction after 
three wave periods (i.e, t = 3 x ^ ), for ^ad — 80, with a fully refined grid using 32 cells. 
The solid line represents the analytical solution. The agreement between the numerical and 
the analytical solution is excellent (see the amplitude of the error on the figure), even after 
the wave amplitude has decreased by a factor of about 2. 

In order to check for the numerical diffusion as explained in Equation (80) we need 
to perform the same simulation using less cells for the numerical diffusivity (r) num ) to be 
not negligible compared to the physical diffusivity (t)ad)- The profiles of B lx , V\ nxi B ly and 
Vi n y along the z direction after five wave periods (i.e, t = 5 x ^ ) for 7^ = 80 with a 
grid of 16 cells is represented Figure 12. The solid lines represent the analytical solutions 
either without taking into account the numerical diffusivity (the not corrected curves), 
or correcting the damping factor according to Equation (80) (the corrected curves). The 
agreement between the numerical and the analytical solution taking into account numerical 
diffusivity is excellent (see the amplitude of the error on the figure). 

The standing Alfven waves test We now start the simulation from an initial 
perturbed state obtained by adding two propagating waves in opposite directions with the 
same damping s r , and let the system evolve. Figure 13 displays a snapshot of the evolution 
of B lx , Vi nx , B^, Vi n y, P an d P along x, after about 4 periods (in order for vi nx and V\ ny 
to be greater than zero) for 7^ = 80. The excellent agreement between the numerical and 
the analytical solution is confirmed. 
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As previously, we determine the time evolution of the pressure thanks to Equation (82) 

2s r t ( s r cos(2sit) + Sjsin(2sjt)\ s r 



P = Pintt + (7 - l)k'r]AD 



+ e 



|s| 2 J 



(84) 



Following Choi et al. (2009) it is interesting to study the time variation of the magnetic 
field in the z direction, B x , as represented Figure 14. The analytical solution is represented 
by the solid line while the dotted line represents the error and the squares the simulation. 



2.4-4- Convergence order 

We tested the evolution of the precision of the implementation of ambipolar diffusion 
by examining the evolution of the error with the level of refinement, i.e with the mesh size 
Ax, for Alfven standing waves and the Barenblatt test. The error e is defined here as the 
maximum difference between the analytical values and the numerical solution, corrected by 
the damping factor for Alfven waves. The error against the cell size follows a power-law, at 
least in the range studied here (up to 10 periods of the wave). For the standing waves we 
find 

e oc Ax 3 . (85) 

For the Alfven propagating waves 

e oc Ax 2 . (86) 

For the Barenblatt test 



e oc Ax 2 . 



(87) 
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A log-log plot of the error as a function of cell size Ax for different times is shown 
on Figure 10 for Alfven standing waves and propagating waves, and on Figure 7 for the 
Barenblatt test. Note that the evolution of the error follows the power laws found through 
the modified equation study, in Section 2.4.3. 

2.4-5. Estimate of the numerical drift coefficient of ambipolar diffusion 

As seen in § 2.4.3, the dissipation of Alfven waves is slightly larger than expected 
according to the analytical values. The spurious dissipation due to the numerical scheme 
can be estimated as: 



I , yKj^j 

^Imes ^lAD 'fnum 

where 7 mes is the value measured in the numerical simulation, with 7 mes _1 = — fr^r, and 
Inum is the drift contribution due to numerical dissipation. Another way to proceed is 
to set 'j ad = oo, to examine how the Alfven waves dissipate, and then to estimate ■y num 
1 — — fl^r- Both methods give about the same value for j nurn . For a level of 
AMR refinement of 2 4 , we get 7„ um _1 = 3 x 10~ 3 ; for 2 5 , 'y num ^ 1 = 5 x 10" 4 and for 
2 6 , Tn.um -1 = 6 x 1CT 5 , to be compared with 7ad -1 = 0.0125 or 0.033 for the present 
simulations. As expected, the better the resolution, the smaller the numerical diffusion. 

Figure 15 is a plot of the dissipation of Alfven waves with ^ad = oo, as explained 
previously. The red solid line corresponds to the analytical solution corrected with our 
estimate of the magnitude of the numerical diffusion, as explained in Equation (80), while 
the black solid line corresponds to the uncorrected analytical solution (no diffusion). 
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3. Ohmic diffusion 
3.1. Equations 

We now turn to the case of Ohmic diffusion in the MHD equations. Equations (1), (2), 
(5) and (8) remain the same. The energy equation is now: 

dE tot 



+ V. (v(E tot + P tot ) - B(v.B) -E ft xB)=0, (89) 



dt 

where E tot and P to t denote the total energy and pressure: 



E tot = pe + \pv 2 + l -B 2 (90) 



2' 2 
1 

2' 



P tot = ( 1 -l)pe+)-B 2 . (91) 



The time evolution of B reads: 



dt 

The Lorentz force and the Ohmic diffusivity EMF read: 



— = V x (v x B -tjqV x B). (92) 



F Lorentz = (V x B) x B (93) 
E n = -t?qV x B, (94) 



where t]q denotes the Ohmic diffusivity. 



3.2. Computation of Ohmic diffusivity 

Various authors (Machida et al. (2006), Machida et al. (2007), Machida et al. (2008), 
Machida et al. (2009)) have studied the influence of Ohmic diffusion, in particular in 
the context of molecular cloud's collapse. Their work assumes that the heating from 
Ohmic resistivity is negligible, and that the approximation V x (— 77c V x B) ~ tjqAB 
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is valid. We choose a more general framework and do not assume either of these two 
assumptions. We implement in RAMSES non- isothermal Ohmic diffusivity, with the exact 
EMF E n = -rjnV x B 

To compute the term of Ohmic diffusivity we proceed exactly as in § 2.2.1. 



3.2.1. The Ohmic diffusion EMF 
The EMF in the z direction En • e z = — rjn(V x B)^ is to be computed at x { _ i , y •_ i , z k . 

2 J 2 

Since the EMF writes: 

E *»-hj-h* - ~H Ax A~y J' (%) 

it is naturally defined at the right position using the natural definition of the Ohmic field 
components (see Figure 1). r]n is computed at x^i , y,-_i, Zk using the procedure described 

2 "* 2 

in § 2.2.1 to compute 7ad, P and pi. 



3.2.2. The Ohmic diffusion energy flux 

This flux writes Tq = r/n(J x B). As explained in § 2.2.4 the flux has to be evaluated 
on each face of the cell, that is at locations (x i± i,yj, z k ), (x ir y j± i, z k ) and (xi, i/j, z k ±i). 

2 J 2 2 

The computation of J and B at these locations is already explained in § 2.2.4. 



3.2.3. Computation of the time step in presence of Ohmic diffusion 



The characteristic Ohmic diffusivity time step, to, is computed according to 

Aa; 2 

t n = 0.1 x , (96) 

m 
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where, as for the ambipolar diffusion case, the coefficient O.f yields a small enough time step 
to ensure good code convergence. The computational time step is the minimum between £q 
and the time step obtained for the ideal MHD case. 

3.3. Tests for the Ohmic diffusion 

3.3.1. Test of Ohmic diffusivity alone 

We first examine the accuracy of the treatment of Ohmic diffusivity alone. We take 
exactly the same conditions as in § 2.4.1 for ambipolar diffusion. We further assume that 
?7n is constant. In that case the induction equation reduces to a diffusion equation with a 
constant diffusion coefficient, using the divergence-free condition V • B = 0: 

f)T2 

-Of = VnAB. (97) 

The solution to this equation for an initial state given by a Dirac pulse is the well 
known heat diffusion equation which yields a gaussian distribution with a width spreading 
as a oc y/i. This can easily be studied either for a one dimensional pulse (e.g. B y (x), 
B x — 0, B z — 0) or a two dimensional pulse (e.g. B y (x,z), B x — 0, B z — 0). The results 
(setting t]q = 1) are displayed on Figure 16. The agreement between the numerical and the 
analytical results is excellent, always better than about 0.5%. We checked that the results 
obtained on an AMR grid are as good as the ones obtained on a regular grid corresponding 
to the highest level of refinement, and that exactly the same results are obtained for any 
orientation of the magnetic field. 

The evolution of the error as a function of the resolution is represented Figure 17. For 
this particular test (heating equation) the spatial scheme is of order 2: e oc Ax 2 . 

In this case, due to the smoothly varying magnetic field, using a refinement strategy 
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based on V-B is not very efficient: the AMR runs are using typically the same number of 
cells as a fully refined grid (at least for our tests between level 5 and 9). 



3.3.2. C-shock 

Proceeding as in § 2.4.2 we have tested the accuracy of our treatment of Ohmic 
diffusivity for the case of an oblique C-shock. For a stationary shock in the x direction (all 
quantities are supposed to only depend on x) the equations of mass, momentum, energy, 
magnetic field conservation and the condition V.B = read: 

d x (pv x ) = (98) 
d x (pv 2 x + P gaz + ^B 2 y ) = (99) 

d x (pv x v y - B x B y ) = (100) 

d x {{E tot + P tot )v x - (v ■ B)^ - rfrB v d x By) = (101) 

d x (v x B y - v y B x - rj n d x B y ) = (102) 

d x (B x ) = 0. (103) 

This set of equations is solved numerically and provides the benchmark to which the 
simulation with the RAMSES code will be compared to assess the accuracy of the numerical 
treatment in the code. 

We start from a steep function as initial state for the different variables whose values 
are the ones taken at infinity ahead of and behind the shock, respectively, in the frame 
of the shock. These values are displayed in Table 2. For this test the Ohmic diffusivity 
coefficient is set to rja = 0.1. The results are portrayed on Figure 18. As seen in the 
figure, after a transitory regime the shock becomes stationary, as expected. A very small 
drift velocity of the shock front persists, of the order of 0.25% of the minimum value of v x . 
Identical results are obtained for any orientation of the magnetic field and of the initial 
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velocity. 

Such an agreement between the numerical and the analytical solutions, within about 
0.2% (except a few points where it can reach 1%) can be considered as very satisfactory and 
asseses the validity of our treatment when hydrodynamics and Ohmic diffusion are coupled. 

The grid is refined if the gradient of magnetic field, pressure, density or velocity is 
greater than 0.1 (this insures for this test that the error on the AMR grid and on the 
regular grid are about the same). 

3.3.3. Alfven waves 

Proceeding as for the ambipolar diffusion study we have examined the behaviour of 
propagating Alfven waves as well as of standing waves in an non-isothermal ionized plasma 
in the case of Ohmic diffusion. Lesaffre and Balbus (2007) derived analytical solutions for 
the general case of MHD flows with shear, non-zero resistivity r]Q, viscosity and cooling. In 
the absence of shear and rotation, these authors showed that torsional Afven waves are a 
solution for such flows. 

Following closely the notations of Lesaffre and Balbus (2007), the unperturbed state in 



Variable 


P 


v x 


Vy 


B x 


By 


P 


Pre-shock value 


0.4 


3 





2 


2 


0.4 


Post-shock value 


0.71084 


1.68814 


0.4299 


V2 
2 


1.43667 


1.19222 



Table 2: Initial state used to generate an oblique C-shock, as described in § 3.3.2. 
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both studies for a wave propagating in the x direction is defined as: 

P = 0.625, pan = 1, Poi = 1, 

v te = o, v^ = o, v 0z = o, 

-Box = 0, = 0, B 0z = 1. 

For the propagating wave study, the perturbed state is chosen such as 5b x = 1 and 
Sb y = iSb x (we necessarily have 5b z = 0). Furthermore, we have 5p = (constant density), 
but the pressure varies with time, so that SP ^ (see Lesaffre & Balbus 2007). In the 
absence of shear, viscosity and rotation, the relation between the perturbed magnetic field, 
Sh = (5b x yi + 5b y y)e st+ ' lkz , and the perturbed velocity, 5u, reads 

s5u = i—5b, (104) 
P 

with s the wave angular frequency and k the wave number. 

The time evolution of the gaz pressure P is governed by the equation: 

dt(-^-r) + V ■ (-^-H = -PV ■ (<5u) + m J 2 , (105) 
7 — 1 7 — 1 

with 7 the adiabatic coefficient of the gaz and J = V x B the current. 

Since Su only depends on z and has components only in the x and y direction, 
div5u = 0. We finally get 



dtP= ( 7 -l)^ J 2 - (106) 
The solutions of the dispersion relation read: 

* = -\ ± ^/(f) a -*M, (107) 

with r) = /c 2 77n- A value 77^ = 5 x 10 -3 yields a moderate damping, with s r = —9.8696 x 10 -2 
and Si = ±1.9844, whereas a value rj n = 2 x 10~ 2 produces a stronger damping, with 
s r = -3.9478 x 10" 1 and s { = ±1.9473. 
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Estimating numerical diffusion Proceeding exactly as in § 2.4.3 we can derive the 
leading order error term in the ideal MHD scheme for Alfven standing, and propagating 
waves. 



Propagating waves We start the simulation from an initial perturbed state with 
5b = 5b x .TZe (e ikx (x + iy)) and Su = ^Ke ( V fcx (x + fy)). The time evolution of the 
pressure is: 

P=P M+ (7 ~ i y 2 *V -'-l). (108) 

For the propagating wave test, we have arbitrarily chosen > and let the system 
evolve from the initial state. Figure 19 portrays a snapshot of the evolution of 8b x , 
5u x , Shy, 8u y , p and P along the ^-direction after five wave periods (i.e, t = ^ 2L ), for 
rjn = 5.10" 3 . Once again, the agreement between the numerical and the analytical solution 
is very satisfactory, at most of the order of a few percents. 



Standing waves As for the ambipolar diffusion, we start the simulation from an initial 
perturbed state obtained by adding two propagating waves with opposite values of and 
the same value of s r , and let the system evolve. The evolution of the pressure is (in real 
notation): 



P =/>,„,« + (7 - + ^iSiM^M) - iy . (109) 



Figure 20 shows a snapshot of the evolution of 5h x , 5u x , 5h y , 5u y , p and P along z after 
three wave periods (t = 3x2,7r ), for t)q = 5 x 10~ 3 . As seen, once again, the agreement 
between the numerical and the analytical solution is very good, of the order of or better 
than a few percents. 



3.3.4- Convergence order 



We tested the evolution of the precision of the implementation of Ohmic diffusion by 
examining the evolution of the error with the level of refinement, i.e with the mesh size 
Ax for Alfven standing waves and the Barenblatt test. The error e is defined here as the 
maximum difference between the analytical values and the numerical solution, corrected 
by the damping factor for Alfven waves, and the error at the center of the box for the 
Barenblatt test. The error as function of cell size follows a power-law, at least in the range 
studied here (up to 10 periods of the wave). For the standing waves we find: 

eccAx 3 . (110) 

For the Alfven propagating waves 

eccAx 2 . (Ill) 

For the Barenblatt test 

eccAx 2 . (112) 

A log-log plot of the error as a function of cell size Ax for different times is given on 
Figure 22 for Alfen standing waves, and on Figure 17 for the barenblatt test. The behavior 
of the error for the propagating waves differs if the mesh size is coarse or fairly refined ( 
e oc Ax 1 or e a Ax 2 respectively). For the Barenblatt test, the error scales as ~ Ax 2 . In 
the case of Ohmic diffusion, Equation (97) reduces exactly to the Heat equation, whereas 
in the case of ambipolar diffusion, Equation (54) reduces to a non-linear diffusion equation. 
The error in the two cases scales as ~ Ax 2 . 
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3.3.5. Estimate of the numerical drift coefficient of Ohmic diffusion 

As seen in Section 3.3.3, the dissipation of Alfven waves is slightly larger than expected 
according to the analytical values. The spurious dissipation due to the numerical scheme 
can be estimated as: 

Vmes = Vn + Vnum, (113) 

where r\ mes is the value measured in the numerical simulation, with r\ mes = — 2Sr <™ m ; and 
Vnum is the drift contribution due to numerical dissipation. Another way to proceed is 
to set tiq = 0, to examine how the Alfven waves dissipate, and then to estimate r] num as 
Vnum — — 2Sr ^ um m Both methods give about the same value for r] num . For a level of AMR 
refinement of 4, we get r\ num — 1. x 1CT 3 ; for 5, r) num — 1.x 10~ 4 and for 6, r] num — 1.1 X 10~ 5 , 
to be compared with rj^ = 0.005 or 0.02 for the present simulations. As expected, the better 
the resolution, the smaller the numerical diffusion. 

4. Conclusion 

In this paper we have described a numerical method to implement the treatment of the 
two important terms of non-ideal MHD, namely ambipolar diffusion and Ohmic dissipation, 
into the multi-dimensional AMR code RAMSES. For ambipolar diffusion, we have used a 
single fluid approach, which is valid when the Lorentz force and the neutral-ion drag force 
are comparable, corresponding to a domain of strong coupling between the fluid and the 
field lines. The situations where such an approximation can be made are numerous, of 
which cloud collapse or certain protoplanetary disks are two typical examples. The accuracy 
of our numerical resolution of the MHD equations was examined by performing a diversity 
of tests, for which either analytical or benchmark solutions exist. For both ambipolar and 
Ohmic diffusion, we first explored the purely magnetic case, without any coupling to the 
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hydrodynamics. For ambipolar diffusion, this was done by comparing the evolution of a 
Dirac pulse to the solution provided by Barenblatt while for Ohmic diffusion, the solution 
is confronted to the well known heat diffusion equation. In a second step, we studied the 
full MHD case (coupling the fluid to the magnetic field) by considering first an oblique 
shock, and then the behavior of propagating and standing Alfven waves. For all these 
tests the solutions obtained with our method show excellent agreement with the analytical 
predictions, typically within a few tenths of a percent on average, showcasing the validity 
and the robustness of our method. We have also carefully analyzed the main source of 
numerical error using the Modified Equation framework. In order to estimate the spatial 
resolution that is required to model non-ideal MHD effects reliably. This opens the avenue 
to a vast domain of astrophysical applications, in particular cloud collapse, pre-stellar 
core formation and protostellar disks where ambipolar and Ohmic diffusion processes are 
believed to play a dominant role. Such astrophysical applications of the non-ideal MHD 
equations with RAMSES will be explored in forthcoming papers. 
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A. The Barenblatt-Pattle solution 



Following Grundy & McLaughlin (1982), the solution of Equation (54) in general form 



> dBy 

- at 



V .(B^S/By)), where (3 depends on the problem, is: 



At a 



IB- 1 



if r < rj t 5 



(Al) 



if r > rj t s 

With [i the dimensionality of the problem, the various constants are defined as follow: 

-fJL 



a 



5 



2 + nf3 
1 

2 + /J/3 



,c2 



B y0 (x)dx = Vo i+2/(3 (-S(3f 



1 r(| / i)r(i//3 + i) 



(A2) 
(A3) 

(A4) 
(A5) 



B. Semi-analytical solution for the isothermal C-shock 



Following Mac Low et al. (1995), in the isothermal case with a constant ion density, we 
reduce the set of MHD equations to: 



pv 2 x + P + § =<h 

pV X Vy ~ ByB X C, 

b 2 - b 2 =2A 2 (D - l)(D~ l - M~ 
(IT 2 - M~ 2 )L^- =-(b 2 + cos^ 1 



(Bl) 
(B2) 
(B3) 

(B4) 



6 _£>(**°cos0 2 + siii0) 

with Ci and Co derived from the initial state, A = — the Alfven Mach number, and M = - 
the Mach number; 9 = 45° is the angle between the magnetic field and the velocity field; 
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and D = and 6 = ^ are the dimensionless density and magnetic field. 



C. Semi-analytical solution for the non-isothermal C-shock 

Following Duffin & Pudritz (2008) and Wardle (1991), and reminding that the set of 
equations is not exactly the same as ours, we solve the set of equations: 

db >jADpioA 2 r 



dx v*b 



(CI) 



l-^r n p\dp lADpior { I 7 S n + sin9\ 

-V 1 (C2) 



(7 - l)r n J dx v s \r n ^ 
b-b 2 



S n = -^coe J 9 (C3) 



1 



(P - Po) - 



V 2A 2 > 

b 2 + cos 2 6 



br n (S n + sin 9) + cos 2 9 



(C4) 
(C5) 



r = 1 - — (C6) 

P B 

where the dimensionless quantities are p = - n " 2 , b = the velocities w nx = 



fnj/ = s "g° Vs = Po and feo are the initial dimensionless pressure and magnetic fields, 
= 45° is the angle between the pre-shock velocity and the magnetic field, and A — ^f- the 
Alfven Mach number. 
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Fig. 10.— Evolution of the error e = l ^ numcrical ^ analyticair with the mesh size Ax 



<-?V (-Bynumcrical— -Bj/analytical ) 2 

n=l N 

for Alfven standing waves (left plot) and Alfven propagating waves (right plot) at different 
times. The dashed lines correspond to two slopes: e oc Ax 3 for the standing waves and 
e oc Ax 2 for the propagating waves. 
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Fig. 11. — The propagating Alfven waves test with ambipolar diffusion {^ad = 80) after 
about five periods. The simulation is represented by squares, while the solid-line is the 
analytical solution. The dotted line is the relative error. We use for this test a fully refined 
Cartesian grid with 32 cells. 
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Position Position 



Fig. 12. — The propagating Alfven waves test with ambipolar diffusion {^ad = 80) after 
about five periods. The simulation is represented by squares, while the solid-lines are the 
two exact solutions (taking into account or not the effect of numerical diffusion according to 
Equation (80)). We use for this test a fully refined Cartesian grid with 16 cells. 
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Fig. 13. — The standing Alfven wave test with ambipolar diffusion ('Jad = 80) after about 
four periods. The simulation is represented by squares, while the solid-lines are the exact 
solutions. The dotted lines represent the relative error. We use for this test a fully refined 
Cartesian grid with 32 cells. 
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Fig. 14. — Time evolution of v < i?x 2 >, the root-mean-square of the magnetic field in the 
x direction at the center of the box, for Alfven standing waves with ambipolar diffusion 
(lAD = 30). The squares are the result of the simulation and the solid line is the analytical 
solution. 
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Fig. 15. — Plot of the magnetic field without ambipolar diffusion: '-/ad = oo. The black 
solid line shows the analytical solution of the unperturbed Alfven wave, while the red solid 
line shows the analytical solution with numerical diffusion taken into account (corrected as 
explained in Equation (80) for a level of refinement of 4). 
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Fig. 16. — Test for Ohmic diffusion only, assuming a Laplacian. The upper panels are 
snapshots for the ID test, with an AMR grid with levels from 5 to 7, at times t = 1.10 -3 on 
the top left and t = 1.10~ 2 one the top right panel. The solid lines are the analytical solution, 
while the dashed lines are the relative error. The lower panels represent the 2D test, on a 
fully refined grid up to level 5, with a contour snapshot on the right and a transverse cut on 
the left, at t = 5.10 -3 : the symmetry is well preserved. 
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Fig. 18. — Non-isothermal oblique shock with Ohmic diffusion (77^ = 0.1). Same caption 
as in the previous figures. The level of refinement is from 5 to 7. 
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Fig. 19. — Alfven propagating waves after five 
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Fig. 20. — Alfven standing waves after four periods and a half, for r) m< i = 5 x 10 3 . 
caption as in the previous figures. 
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Fig. 21. — Time evolution (expressed in units of the period, ^ 2E ) of v < Bx 2 >, the root- 
mean-square of the magnetic field in the x direction, for Alfven standing waves. In this case, 
r]n = 2x 10- 2 . 
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Fig. 22.— Evolution of the error e = l ^ numcrical ^ analyticair with the mesh size Ax 
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for Alfven standing waves (left plot) and Alfven propagating waves (right plot) at different 
times. The dashed lines correspond to two slopes: e oc Ax 3 for the standing waves and 
e oc Ax 2 for the propagating waves. 
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